Nonlinear elastic model for faceting of vesicles with soft grain boundaries 
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We use an elastic model to explore faceting of solid-wall vesicles with elastic heterogeneities. We show that 
faceting occurs in regions where the vesicle wall is softer, such as areas of reduced wall thicknesses or concen- 
trated in crystalline defects. The elastic heterogeneities are modeled as a second component with reduced elastic 
parameters. Using simulated annealing Monte Carlo simulations we obtain the vesicle shape by optimizing the 
distributions of facets and boundaries. Our model allows us to reduce the effects of the residual stress generated 
by ciystalline defects, and reveals a robust faceting mechanism into polyhedra other than the icosahedron. 



Nature uses various vesicular structures to compartmental- 
ize and differentiate matter in order to achieve a desired set 
of functions. Depending on their origin and function, vesicles 
can take a wide variety of shapes and sizes with very differ- 
ent physical and chemical properties. A prototypical example 
of a faceted structure is the capsid of large viruses, for which 
an icosahedral shape was elegantly described by Caspar and 
Klug Ol]. Next to these regular polyhedra, bacterial micro- 
components such as carboxysomes [2] (found in cyanobacte- 
ria (so, for example), are known to take on irregular faceted 
shapes. X-ray diffraction of small (~ 30nm diameter) single- 
walled DPPC vesicles found |4?] that below the thermal tran- 
sition they facet into polyhedra. Thick-wall nested fullerenes 
are also known to sometimes take faceted shapes |5]. Faceted 
shapes were recently observed [6i by coassembling oppositely 
charged amphiphiles. It has been argued that under the right 
conditions electrostatic correlations can lead to faceting fTt] or 
even a complete collapse of the structure [8]. 

The existence of facets separated by sharp edges that have 
been observed in experiments can be understood in terms of 
the vesicle wall being, at least in pait, solid [4] or as a result 
of excess amphiphiles exhibiting spontaneous curvature being 
present in sufficient quantity |9nl^. Molecules within a solid 
wall have fixed connectivity and the vesicle can sustain shear, 
bending, and stretching deformations. The microscopic struc- 
ture and properties of the ordered phase is intimately related 
to the vesicle topology. Very interesting phenomena occur if 
the local molecular order favored by interactions is not com- 
patible with the global topology and thus cannot be realized 
simultaneously throughout the whole surface [11,]. For exam- 
ple, the preferred planar ordeiing of identical paiticles inter- 
acting with a spherically symmetric pair potential is a hexag- 
onal lattice, where each site has exactly six equidistant near- 
est neighbors. In other words, it is possible to tile a plane 
with equilateral triangles. However, this is not the case for a 
sphere, and topological defects (i.e., sites with coordination 
z ^ 6) are necessarily present. For a triangular tessellation 
of a sphere, Euler's theorem ensures that there are at least 12 
sites that are five-coordinated [11]. These disclination defects 
repel each other and take a conformation that maximizes their 
mutual separation, thus positioning themselves at the vertices 
of an inscribed icosahedron. 



The stretching en ergy of a five-fold disclination on a plate 
of radius r is cx lll2ll . If the plate is allowed to buckle out 
of plane it can reduce the stretching energy at the expense of a 
bending penalty that is oc log r [ 12] and form a cone with the 
apex centered at the defect. On a sphere these two energies 
follow the same scaling laws, albeit with different prefactors 
[HU. If the sphere is sufficiently large it can lower its energy 
by buckling into an icosahedron via a similar conical defor- 
mation seeded at the twelve defects [14]. Recently we have 
shown that if a vesicle is assembled of multiple molecular 
species with different elastic properties, depending on the rel- 
ative ratio and the line tension of the components, it can take a 
number of regular and irregular polyhedral [15| or Janus-like 
shapes [ 16]. 

These results suggest that a more general buckling mech- 
anism is possible when the vesicles are not elastically ho- 
mogeneous. However, in a previous work fisl we studied 
two-component elastic shells with the twelve five-fold defects 
fixed to the corners of an inscribed icosahedron. Here we 
show that faceting occurs even if such restrictions are lifted. 
In this paper we introduce a nonlinear elastic model to study 
faceting of a vesicle with solid domains connected with soft 
boundaries. We argue that material heterogeneities, such as 
local variations in the wall thickness or due to defects in the 
crystalline order, can lead to local softening of the vesicle wall 
and provide a pathway for lowering the energy. In order to sat- 
isfy both, the preferred local flatness of the solid domains and 
the curved global geometry the curvature is focused along the 
softer domain boundaiies, which leads to faceting of the entire 
vesicle. 

We assume that the vesicle wall is thin compared to its ra- 
dius, R, and can be represented as its two-dimensional mid- 
surface. Such an approximation is indeed justified as the vesi- 
cle wall is typically a few nanometers thick compared to the 
radius that ranges from tens of nanometers to microns. The 
midsurface can be parametrized with two parameters and 
s^, as a locus of points r = r (s^, s^). The metric tensor is 
defined as — daV ■ 9^r, where a, (3 = 1,2. The second 
fundamental form, related to the surface curvature, is given 
by bai3 = dapY ■ n, with n being the unit normal. The elastic 
energy of the midsurface can be written as ifpTl 
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gajs) is the strain tensor, g^j^ is a ref- 
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determinant, h Ris the vesicle thickness and the summa- 
tion over pairs of repeated indices is assumed. Compared to 
a recent study |fl8|] . here we follow Koiter's arguments lfl9ll 
and retain only terms consistent with the Kirchhoff-Love as- 
sumptions [201 of a negligible normal stress. The rank-four 
elastic tensor A"^^^ is determined by the material proper- 
ties. If the material is isotropic A"^''^ depends only the 
Young's modulus Y and Poisson's ratio and A"^''^ — 

^ {j^r^r' + r^t') , with r'^g^^ = -5^ and 61 be- 
ing the Kronecker delta symbol. 

The first term in Eq. ([T) represents energy associated with 
stretching deformations while the second term accounts for 
bending. For an isotropic material the stretching and bending 
energies simplify to 

= (3) 

with = 'g'^^Ua-y and = 'g^'^bay With the mean 
curvature H = ^b" and the Gaussian curvature K = 
det(b^), the bending energy can be rewritten as Eb = 
J dAK (2i/2 - (1 - I/) iC), where k = ^ ^j^^ is the bend- 
ing rigidity. The prefactor of the Gaussian curvature term 
is negative and proportional to k, thus suppressing saddle- 
like (K < 0) conformations. For a homogeneous vesicle 
with fixed topology the Gauss-Bonnet theorem ensures that 
J dAK — const, and the Gaussian curvature term can be 
omitted. However, for a multicomponent vesicle, as is the 
case here, this term must be retained. 

Variation of Eq. ([T]i with respect to u^p and bap leads to 
a set of nonlinear partial differential equations for the me- 
chanical equilibrium [Is', 120^]. Instead of directly solving 
those equations, which is a formidable task even in simple 
geometries, we opt for a numerical minimization of the en- 
ergy, Eq. ([Til, via Monte Carlo (MC) simulations. The vesicle 
is represented as a discrete triangular mesh with two types 
of triangles, soft and hard. From Eqs. (|2]i and ([Sj it is evi- 
dent that for a fixed Young's modulus Y and Poisson's ratio v 
both bending and stretching moduli depend only on the vesi- 
cle wall thickness, h. Therefore, for simplicity, we assume 
that the elastic properties of a component are determined by 
its thickness, with the soft component being thinner This al- 
lows us to explore the phase behavior as a function of only 
two parameters, relative thickness 77 — j^''^ and the fraction 
/ of the soft component. We keep in mind, however, that in 
the experimental systems the elastic properties are likely to be 
determined by the local organization of the molecules rather 
than by the thickness variations. The discrete stretching en- 

ergy is M Ef^- = (t^ (Tr^)' + TrF^^, 

where F ^'g — I is the Cauchy-Green strain tensor. At 
is triangle area and the sum is carried out over all triangles. 
Matrices 'g and g are the discrete versions of the reference 
and actual metric tensors, respectively, whose elements are 
the scalar products of the two vectors spanning each triangle 



before and after the deformation. The bending energy of the 
discrete mesh is [121 kJ^t t ~ ' ^j) with ii = -^k 
lfl2[ I22I1 and the sum runs over aU pairs of the nearest neigh- 
bor triangles. We point out that the last expression is actually a 
discrete version of the continuum energy / dAn {2H^ - K) 
lfl2ll . Therefore, we also need a discrete expression for the 
Kv J dAK term. Gaussian curvature at a vertex i is ll23ll 
Ki = KdAi = 27r — where is the angle of the 

adjacent triangle T at j and Ai — ^ J^t is the associated 
vertex area. Finally, |23| kv f dAK — > kvJ^t X"' 
where the /-sum runs over the three vertices of triangle T. 

The surface mesh was constructed by building a triangu- 
lation with vertices randomly (but evenly) distributed on 
a sphere. The mesh with Nt triangles was generated with 
the 3D Surface Mesh Generation package in the CGAL li- 
brary ||24| . In order to ensure that the results are insensitive to 
the discretization details we have performed independent runs 
starting from different initial configurations with 2 x lO'^ to 
5 X 10'^ vertices and 4 x 10'^ to 10"' triangles. Furthermore, 
we have validated our results using regular triangulations of 
comparable sizes constructed according to the prescription 
introduced by Caspar and Klug [1]. = fNt, where 

/ = 0.01 . . . 0.4 is the fraction of randomly chosen triangles 
were designated as soft. For a given random mesh we have 
explored at least two different initial random distributions of 
the soft triangles, removing a possible bias caused by a partic- 
ular choice of the initial distribution of the components. The 
reference metric 5 is an input parameter in our model and was 
chosen to be that of the initial triangulation, which removes 
the instability toward buckling into an icosahedron discussed 
in Ref. |14]. However, some residual stress is still present 
as the stretching energy is minimized in the spherical config- 
uration while the bending energy favors flat faces. In order 
for the bending energy to win, the vesicle wall has to be suf- 
ficiently thick. From Eqs. (|2|i and (|3]l bending and stretching 
energies are comparable at /i « a, where a is the unit length 
set by the average edge length of a triangle. Our model implic- 
itly assumes that each triangle contains a sufficient number of 
microscopic degrees of freedom such that the molecular de- 
tails are of no importance. In an typical amphiphilic system 
a 5nm sets the length scale down to which the contin- 
uum description is applicable. We set hfiard — a (a ^ R 
and the thin plate approximation is valid) and note that in the 
hfiard — > limit the vesicle would remain spherical regardless 
of the thickness ratio 77 — hh<i^_ Therefore, in our model a fi- 

nite thickness is essential for faceting to occur Setting g to be 
flat leads to interesting effects which will be addressed else- 
where. In order to find the optimal vesicle shape and the dis- 
tribution of components simulated annealing MC simulations 
were performed. A MC step consisted of two moves: (i) A 
random displacement of a vertex was attempted, followed by 
(ii) an attempted swap of the component types of a randomly 
chosen pair of triangles. Both moves were accepted accord- 
ing to the Metropolis algorithm. The system was heated up 
and cooled down using a linear cooling protocol followed by 
three consecutive exponential cooling runs. For each run, the 
configuration found to have the lowest energy was recorded. 
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Figure 1. (color online) Snapshots of typical vesicle shapes as a func- 
tion of the thickness ratio rj = hhard/hsoft ~ 1-5, 2.0, 2.5 and the 
fraction of the soft component / = 0.06, 0.19, 0.24 and 0.32. The 
soft component is shown in blue (dark) while the hard component is 
white (light) and slightly transparent. For clarity, width of the soft- 
component region has been enhanced. 



A typical run consisted of 4 x 10^ MC sweeps. Note that the 
triangle type swap move is just a tool optimize the component 
distribution and does not correspond to the actual redistribu- 
tion of the molecules that occurs during the vesicle assembly. 
The actual kinetics of the self-assembly process is complex 
and beyond the scope of this work. 

We set Poisson's to = i, measure energies in units of 
Yh/a^ and measure mean and Gaussian curvatures in units of 
1/a and 1 /a^, respectively. Note that the energy landscape is 
complex and the simulated annealing is unlikely to find actual 
ground state configurations. However, all obtained structures 
are qualitatively reproducible and are thus referred to as typi- 
cal. In Fig. [T]we present snapshots of typical vesicle shapes 
for a range of fractions / of the soft component and relative 
thicknesses, 77 — ^^'"^'^ . For / < 0.1, the soft component 
forms elongated, mutually disconnected ridges on the vesicle 
surface. The vesicle remains nearly spherical with slight dis- 
tortions near the ridges. As / increases to w 0.2, ridges begin 
to merge and facets develop. The onset of faceting is not sharp 
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Figure 2. (color online) Distribution of stretching Es (left) and bend- 
ing Eb (right) energies for rj = 2.0 and / — 0.25. Black lines indi- 
cate the soft boundaries between facets. For visualization purposes, 
a lattice dual to the actual triangulation is shown. 



and appears to be sensitive to rj. If the amount of the soft com- 
ponent is further increased the number of facets increases and 
they become smaller in size. We speculate that for / > 0.5 
this would no longer be the case and the soft component would 
form islands; however, we have not explored this regime as it 
is not applicable to the experimental systems of present inter- 
est. The hard-component facets are, however, not perfectly 
flat and their curvature depends on 77. For ry < 1.5 the vesicle 
appears almost spherical, while for rj > 2.5 its shape is clearly 
polyhedral. 

In Fig. |2] we show the distribution of the stretching and 
bending energies on the vesicle surface for rj = 2.0 and 
/ = 0.25. As expected, the stretching energy is highest 
along the boundaries where the material is soft. Note that 
the stretching energy is paiticulaiiy large near the kinks in 
the soft boundary suggesting that such structures are probably 
not energetically prefeiTed. Relaxing these kinks requires a 
coherent global redistribution of the components that is very 
hard to achieve in most numerical optimization schemes. 

In Fig. [3] we show the distribution of the mean curvature 
for a vesicle with / = 0.25 and 77 = 2.0 with two distinct 
peaks. The sharper peak at low H corresponds to the facets, 
while the weaker broader peak at a larger H represents the 
curvature at boundaries. A projection plot of H clearly shows 
that the high curvature is condensed along lines. Distribution 
of the Gaussian curvature is shown in Fig. |4]with a sharp peak 
at zero. Although if = on the facets, H > indicating 
that the facets are not flat but locally cylindrical. Gaussian 
curvature is localized along the soft grain boundaries as can 
be seen in the projection plot in the inset in Fig. H] We note 
that a number of vertices have K < 0. Such regions locally 
resemble saddles and are located along the soft-component 
ridges. A detailed analysis of this effect would require a finer 
mesh and a more accurate minimization technique than used 
here. 

We have implemented a general non-linear elastic model 
to analyze faceting of solid-wall vesicles with soft domain 
boundaries into polyhedra other than icosahedra. By choos- 
ing a suitable reference metric state, we have removed the in- 
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Stability toward buckling into the icosahedron seeded at the 
twelve five-fold defects ifTill . While topological defects are 
still present as required by the spherical topology their ef- 
fects are treated implicitly by assuming that they collect at the 
boundary between facets making the vesicle wall locally soft. 
Consequently the long-range effect of the stress produced by 
the defects are substantially suppressed resulting in faceting 
into irregular polyhedra. We believe that our model quali- 
tatively explains the faceting observed in small DPPC vesi- 
cles |4j] and vesicles coassembled by oppositely charged am- 
phiphiles We conclude by noting that other nonspheri- 
cal shapes have also been observed in multicomponent liquid 
vesicles 1I23I1 a nd in vesicles assembled of smectic polymer 
molecules l25l|26ll . 



Figure 3. (color online) Distribution of the mean curvature, H, for 
/ = 0.25 and ri = 2.0 with two peaks corresponding to the hard 
facets (lower H) and soft boundaries (higher H). Inset: Projec- 
tion, (^0-plot of H. Darker colors indicate high curvature regions. 
Blurry distortions near the plot edges are an artifact of the interpola- 
tion method. 



0.81 1 1 r 




Figure 4. (color online) Distribution of the Gaussian curvature, K, 
for / = 0.25 and rj — 2.0. Inset: Projection, (f)6-plot of K. Darker 
colors indicate high Gaussian curvature regions. Blurry distortions 
near the plot edges are an artifact of the interpolation method. 
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